
#include "header.h"

int main(int argc, char** argv) {

 struct grid   m[128];           // array of grid structures
 struct grid  *temp, *left;
 struct SSdata S;                // contains self similar data
 struct INPdata input;            //  contains startup parameters  
 //const double pi = 3.14;
 double	      phase[3], phi[3];
 char         Lname[80];         // runname (global)
 double       dt, d, t; 
 int          k, q, kk, j, ac = argc;  
 char         filename[80];
 char**       av = argv;
 FILE*        fh;
 double       sL , p;              // selfsimilar parameters described later 
 double	      ampl = 1.98;	  // Amplification factor that determines when to make refinement. 
 				  // (ampl > 1, the lower the more accurate solution will be)



 /*-- startup  --*/

 strcpy( input.runname , argv[1] );
 startup( m , &input , &ac , &av );


 //----------- Initializing main loop parameters  ------------//

 kk=0;
 t=0;
 k=0;                                               
 p=0;                                    
 q=0;                                
 dt = (input.lambda)*pow(input.dr, 2);               
 d  = (input.A); 

 //----------- Initialize self-similar parameters ------------//

 phase[0]   = atan( (m[0].vgrid[1])/(m[0].ugrid[1]) );
 phase[1]   = atan( (m[0].vgrid[2])/(m[0].ugrid[2]) );
 phase[2]   = atan( (m[0].vgrid[3])/(m[0].ugrid[3]) );
 phi[0]     = sqrt( pow(m[0].ugrid[1] , 2) + pow( m[0].vgrid[1] , 2) );
 phi[1]     = sqrt( pow(m[0].ugrid[2] , 2) + pow( m[0].vgrid[2] , 2) );
 phi[2]	    = sqrt( pow(m[0].ugrid[3] , 2) + pow( m[0].vgrid[3] , 2) );

 S.A        = 0.;
 S.Lt       = 0.;
 S.L        = 2.206/(input.A);
 S.Sstime   = 0.5*dt/(S.L * S.L);
 S.Phtime   = 0.;                      // 0.5*(S.L)*(S.L)*dt;
 S.phase    = phase[1];
 S.phase_rr = (phase[0] - 2.*phase[1] + phase[2])/((m[0].dr)*(m[0].dr));
 S.phi 	    = phi[1]; 	
 S.phi_rr   = (phi[0] - 2.*phi[1] + phi[2])/((m[0].dr)*(m[0].dr));	

 strcpy( Lname , input.runname);
 strcat( Lname , ".out");
 fh       = fopen( Lname , "w");

 fprintf(fh,"#  1.Time\t\t 2.|psi|\t\t 3. |psi|_rr\t\t 4. phase_rr\t\t 5. phase\n\n"); 
 printf("Physical Time \t\t  SS L \t\t Number of Particles \n\n");

 //--------------------- Main loop ---------------------------//  
 while (S.Phtime <= input.tmax){  
   if (k == -1) {	
      non_local(0.,m);
      break;
   }
   //-------------------- File Write --------------------------//
   
   if (S.Phtime+0.5*dt >= p){   
     temp = m;
     while ((temp->lpointer)!=NULL) {
         temp = temp->lpointer;
     }
     sprintf(filename, "%s_%04d.bin", input.runname, q);
     printf("%.10e \t %.10e \t %.10e \n", S.Phtime, S.L, Nparticle(m));
     fclose(fh);
     //----------------------------------------------------------------//
     //      Uncomment to enable writing files with profile                                                            //
     //----------------------------------------------------------------//
                                                                      
     fh=fopen(filename, "wb");                                       
     rhs(0., m);       
     left = temp;
     fprintf(fh,"#\tFile contains all grid data for restarts\n");
     fprintf(fh,"#\t%d\n", input.Ngrids);
     for (j=0; j<input.Ngrids; j++){
       fprintf(fh, "#\t%d\n", temp->Np);
       temp = temp->rpointer;
     }
     temp = left;
     
     fprintf(fh,"#time = \t%f\n",S.Phtime);
     fprintf(fh,"#\t1. Radius\t\t\t\t2. abs(psi) \t\t\t\t3. re(psi)\t\t\t\t4. im(psi)\n");
     fprintf(fh, "\n");

                                             
     while (temp!=NULL) {                                        
       fprintf(fh , "\n");                                        
       dump(fh, temp, &(S.Phtime) );                            
       temp = temp->rpointer;                                       
     }                                                                
     fclose(fh);                                                       
     //-------------------------------------//
     fh=fopen( Lname , "a");
     p += (input.tinc);
     q++;
   }


   //----------------------------------------------------//
   //                                 Runge Kutta 4 step                                        //
   //----------------------------------------------------//
  
   temp = m;
   while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}
   left = temp;   
   RK4(m, &dt, 0.);
   k++;

   //-----------------------------------------------------//
   //                           Recalculate self-similar parameters                               //    
   //-----------------------------------------------------//
   
   phase[0]   = atan( (left->vgrid[1])/(left->ugrid[1]) );
   phase[1]   = atan( (left->vgrid[2])/(left->ugrid[2]) );
   phase[2]   = atan( (left->vgrid[3])/(left->ugrid[3]) );
   phi[0]     = sqrt( pow( left->ugrid[1] , 2) + pow( left->vgrid[1] , 2) );
   phi[1]     = sqrt( pow( left->ugrid[2] , 2) + pow( left->vgrid[2] , 2) );
   phi[2]     = sqrt( pow( left->ugrid[3] , 2) + pow( left->vgrid[3] , 2) );

   sL  = S.L;
   S.L      = 2.206/phi[1] ;
   S.Lt     = (S.L - sL)/dt; 
   S.A      = - S.L * S.Lt;
   S.Sstime = S.Sstime + dt/(S.L * S.L);
   S.Phtime = S.Phtime + dt;
   S.phase    = phase[1];
   S.phase_rr = (phase[0] - 2.*phase[1] + phase[2])/((left->dr)*(left->dr));
   S.phi      = phi[1];
   S.phi_rr   = (phi[0] - 2.*phi[1] + phi[2])/((left->dr)*(left->dr));


   if (k>kk){
     kk+= (input.Noutfile);     
                          
     fprintf(fh , "%.15e\t%.15e\t%.15e\t%.15e\t%.15e\n",S.Phtime, S.phi, S.phi_rr, S.phase_rr, S.phase);
   }

   //---------------------------------------------------------//
   //                       Adaptive Refinement                                                        //
   //---------------------------------------------------------//

   temp = m;
   while ((temp->lpointer)!=NULL) {
         temp = temp->lpointer;
   }
   rhs(0.,m);
   /*if (ampl*d < sqrt( pow( temp->ugrid[2] , 2) + pow( temp->vgrid[2] , 2)   )) {
     d=(sqrt( pow( temp->ugrid[2] , 2) + pow( temp->vgrid[2] , 2)   ));
     addgrid( m, &(input.Ngrids));     
     dt             = 0.25*dt;
     input.tinc     = 0.25*(input.tinc);
     //input.Noutfile = 0.5*input.Noutfile;	
   }
   temp = m;
   while ((temp->lpointer)!=NULL) {
     temp = temp->lpointer;
     }
   */
   //-----------------------------------------------------------//
} 
 printf("\t\tAt the end of physical simulation\n\n" );
 printf("\t\tPhysical Time \t SS Time  \t SS A \t\t SS L \n\n");
 printf("\t\t %e\t %e \t %e \t %e \n\n", S.Phtime, S.Sstime
	    - 0.5*dt/(S.L * S.L), S.A, S.L);

 return 0;
}
